source('00_util_scripts/mod_bulk.R')
library(edgeR)

dss.tespex <- read_delim('~/learn/erv_quant/TEspeX/dsstest/outfile.txt')

lib.tespex <- read_delim('~/learn/erv_quant/TEspeX/dsstest/mapping_stats.txt')

dss.meta <- read_delim('~/append-ssd/alaria2/dss_mice/nxf_input.csv') |>
  mutate(group = str_extract(sample, 'DSS|Water') |> fct_relevel('Water'),
         name = str_extract(fastq_1, 'SRR\\d+'), .keep = 'none') |>
  arrange(name) 

dge.tes <- dss.tespex |>
  mutate(id = str_remove(TE, '#.+'), .before = TE, .keep = 'unused') |>
  pivot_longer(-1) |>
  mutate(name = str_remove(name, '.cleaned')) |>
  pivot_wider(id_cols = 1, names_from = name, values_from = value) |>
  column_to_rownames('id') |>
  DGEList(lib.size = lib.tespex$mapped, group = dss.meta$group,remove.zeros = T)

dge.tes <- dge.tes[filterByExpr(dge.tes), , keep.lib.sizes = T]

dge.tes <- dge.tes |>
  calcNormFactors() |>
  estimateDisp()

dge.tes |> plotMDS(labels = dge.tes$samples$group)

# DEGA -------
## LRT: aggressive -------
lrt.tes <- dge.tes |>
  glmFit() |>
  glmLRT()

lrt.tes |>
  decideTests() |>
  summary()

## QLF: conservative ----
qlf.tes <- dge.tes |>
  glmQLFit() |>
  glmQLFTest()

qlf.tes |>
  decideTests() |>
  summary()

## Treat: conservative test -------
treat.tes <- dge.tes |>
  glmFit() |>
  glmTreat()

treat.tes |>
  decideTests() |>
  summary()

qlftreat.tes <- dge.tes |>
  glmQLFit() |>
  glmTreat()

qlftreat.tes |>
  decideTests() |>
  summary()

# visualization ----------
gene.lrt <- lrt.tes$table |>
  as_tibble(rownames = 'gene') |>
  mutate(padj = p.adjust(PValue, method = 'BH'))

gene.lrt |>
  ggplot(aes(logFC, -log10(padj))) +
  geom_point()

pub.erv <- read_csv('mission/fdx1/publication.erv.csv')

gene.lrt |>
  filter(gene %in% pub.erv$gene) |>
  arrange(padj)

cpm(dge.tes) |>
  as_tibble(rownames = 'gene') |>
  filter(gene %in% pub.erv$gene)
  
